A guided Monte Carlo method for optimization problems 
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We introduce a new Monte Carlo method by incorporating a guided distribution function to the 
conventional Monte Carlo method. In this way, the efficiency of Monte Carlo methods is drastically 
improved. To further speed up the algorithm, we include two more ingredients into the algorithm. 
First, we freeze the sub-patterns that have high probability of appearance during the search for 
optimal solution, resulting in a reduction of the phase space of the problem, a concept inspired by 
the renormalization group equation method in statistical physics. Second, we perform the simulation 
at a temperature which is within the optimal temperature range of the optimization search in our 
algorithm. We use this algorithm to search for the optimal path of the traveling salesman problem 
and the ground state energy of the Anderson-Edwards spin glass model and demonstrate that its 
performance is comparable with more elaborate and heuristic methods. 

PACS numbers: 02.60.Pn, 05.70.Jk, 75.10.Nr, 89.80.+h 



Optimization problems arise in areas of science, en- 
gineering and other fields. Many of them are NP hard 
problems, in the sense that the number of computing 
steps required to solve the problem increases faster than 
any power of the size of the system. The traveling sales- 
man problem (TSP)[1], the spin glass problem[2] and the 
Lennard-Jones microcluster problem[3] are some of the 
examples belonging to this class. Over the years, peo- 
ple have developed heuristic methods that can allow one 
to possibly obtain optimal solutions in some optimiza- 
tion problems, e.g. the TSP. However, these methods 
are usually specially designed for the problem that one is 
seeking to solve. In many cases, these heuristic methods 
will either not be applicable or will take a much longer 
time to locate the optimal solution of the problem. In 
many practical situations, it is more desirable to obtain 
a near-optimal solution quickly than to find the true op- 
timal solution of the problem. Furthermore, it is also 
preferrable to have an algorithm that is simple to use 
and general enough to treat various optimization prob- 
lems. 

Since its introduction about half a century ago, the 
Monte Carlo (MC) method [4] has been widely used to 
treat optimization problems. Despite its simplicity and 
versatility, it is known that Monte Carlo method cannot 
usually give one an acceptable solution since it can easliy 
get trapped in local minima of the problem. In the past 
decades, people have developed new algorithms based on 
the original Monte Carlo method. Simulated Anneal- 
ing(SA) [5] is one of the algorithms that originates from 
the Monte Carlo method. Because of its simplicity and 
versatility, SA has been employed to study varoius sta- 
tistical systems with success. Despite its improvement 
over the original Monte Carlo method, SA also suffers 
from a similar problem, that the algorithm can easliy 
get trapped in local minima. Progress has been made 
to improve SA and Monte Carlo methods further by de- 
velopling more efficient MC sampling algorithms. More 



recent developments include algorithms such as the mul- 
ticanonical method [6] and simuulated tempering [7] . All 
the above mentioned algorithms are meant to give a more 
efficient MC sampling over the original method. 

In this paper, we introduce a new algorithm by incor- 
porating the Monte Carlo method with a guided distri- 
bution function. The advantage of this new algorithm 
is that it is very simple to use and is generally suit- 
able for most optimization problems, whether they in- 
volve continuous or discrete parameter spaces. The es- 
sentials of this algorithm is from both biology and the- 
oretical physics. In the following, we will give the main 
ingredients of this algorithm and apply it to two prob- 
lems, namely, the traveling salesman problem and the 
spin glass problem (SGP). The task of the TSP is to 
find a route through a given set of N cities. The spin 
glass model we are testing here is the three-dimensional 
Edwards-Anderson Ising spin glass model [8] in which 
Ising spins (Si — ±1) are placed on a three dimensional 
lattice and there are nearest neighbor interaction only. 
The corresponding Hamiltonian for this model is 



H 



<i,3> 



JijUiUj 



where in our simulation, the Jy takes on values ± J ran- 
domly and (Ji,(Jj can take values +1 or -1. The task 
here is to find the spin values which minimizes H . As it 
turns out, its performance in these two tests is compa- 
rable to the more elaborate methods and it can obtain 
near-optimal solutions rather quickly. 

A new substance in this algorithm is the addition of a 
guided distribution function to the original Monte Carlo 
method. This introduction of the guided distribution 
function is motivated from biology. Similar to the subject 
of evolutionary programming [9] , this guided distribution 
function is based on the cooperative effort of individu- 
als. To begin with, we have a population of M randomly 



generated solutions. We let each of them evolve indepen- 
dently and according to the conventional Monte Carlo 
rules, i.e. we accept a new solution over its old solution 
when the new solution is better fit than the old solution. 
If a new solution is less fit than the old solution, we pick 
the new solution over the old solution with respect to a 
Boltzmann factor. At this point, it is equivalent to M in- 
dividual runs of the Monte Carlo method. In each of the 
M individual runs, we keep record of its best fit solution 
while the individual solution is evolving. 

After a certain number of Monte Carlo steps, a distri- 
bution function of the M best fit Monte Carlo solutions 
will be recorded. Let us use the TSP as an example. Af- 
ter we perform a preset number of Monte Carlo steps, 
each of the M independent runs has its own best fit so- 
lution, or path. In each of these paths, we record the 
links between any two cities. There are a total of N links 
in each path and a total of MN links of the M best fit 
solutions. One has now a distribution function of the 
probability of appearance among all possible links that 
connects any two cities. This completes our first layer of 
Monte Carlo simulation. In our discussion, a layer of sim- 
ulation always means a set of M individual Monte Carlo 
runs for a preset Monte Carlo steps plus the evaluation 
of the distribution function that we just mentioned. The 
distribution function that we obtain here will be used as 
a guided distribution function in the next layer of Monte 
Carlo simulation of M individual runs. 

In the next layer of Monte Carlo simulation, a set of M 
individual Monte Carlo runs is again performed. An im- 
portant point here is that we also start with M randomly 
generated solutions at the beginning of the simulation in 
this layer. There is however one more criterion here — 
that we pick the links which appear less frequent in the 
guided distribution function and try to change them into 
links with a higher probability of appearance. In practice, 
we pick one of the two cities that are connected with a less 
probable link and connect it to a city where this new link 
appears more often in the guided distribution function. 
The new solution will be compared with the old one to 
decide whether we keep or disgard it according to the con- 
ventional Monte Carlo method. The idea behind this is 
simple. The more (less) often that a certain link appears 
among the M best fit solutions, the more (less) likely that 
it would (not) appear in the optimal solution. One would 
then let a given solution evolve into solutions with links 
of high probability of appearance while avoiding links of 
low probability of appearance. The Monte Carlo method 
still allows enough fluctuation to search for optimal, or 
near-optimal solutions. After a preset number of Monte 
Carlo steps is performed, a new guided distribution func- 
tion will again be obtained for the next layer simulation. 
In principle, the more Monte Carlo steps one performs in 
each layer, and the more layers of simulations one has, 
the better will be the result. 

Two more ingredients can be added into this new al- 
gorithm. The first is the elimination of irrelevant degrees 
of freedom. Again taking the TSP as an example. There 



are a total of — ^ — - possible links among the N cities. 
Some of them have a higher probability of appearing in 
the best fit solutions while others have very low probabil- 
ity. We can elminate some of the links by freezing a link 
between two cities that exceeds a certain probability of 
appearance in the guided distribution function. In prac- 
tice, after the guided distribution function in a layer is 
obtained, one freezes links which have a high probability 
say, 0.8 of appearance before one performs the next layer 
simulation. In this way, one reduces the available degrees 
of freedom in the system in the next layer, thus speeds 
up the search process. This idea is inspired by the renor- 
malization group equation method in statistical physics 
and is similar in spirit to a recent work [10] which has 
proposed to treat optimization problems using the idea 
of renormalization in statistical physics. By freezing high 
probability links, one can obtain near-optimal solutions 
within a much faster CPU time. One can unfreeze the 
frozen links during the simulation to avoid missing the 
optimal solution. 

The last ingredient that we include in this algorithm is 
from an observation in a recent work [11]. It was [11] ob- 
served that the average first passage time to reach the 
global optimum of an optimization problem in Monte 
Carlo simulation was a U-shape temperature dependent 
curve. At the optimal temperature, the avereage first 
passage time to find the global optimal will be the short- 
est. This means that for the optimization problem one 
has, it is best to perform simulation at the optimal tem- 
perature so that one can locate the optimal, or near- 
optimal solution in the fastest CPU time possible. The 
reason for the appearance of such an optimal temperature 
can be understood as follows. At a given temperature, 
the system tends to stay in states near the equilibrium 
value determined by that particular temperature rather 
than going to the global optimum. As the system re- 
laxes toward this equilibrium value, it would fluctuate 
around this equilibrium value. The Gaussian like fluctu- 
ation would then carry the system from its equilibrium 
value to the global optimum as long as there are no insur- 
mountable barriers between them. When the simulation 
is performed at high temperatures, there arc too many 
paths available and it would take the system too long to 
locate the path that would reach the global optimum. As 
the temperature lowers, the number of available paths 
are reduced. When the temperature further decreases, 
the potential barriers and traps make the system more 
difficult to move out of the local optima. There is there- 
fore an optimal temperature at which the two factors 
balance each other and the system approaches the global 
optimum with the shortest time. 

(1) The Traveling Salesman Problem (TSP) 

Since we have already used the TSP as an example 
to explain the implementation of our algorithm above, 
we do not need to repeat the essential steps here but 
instead only state the parameters we used and give some 
comments at the end. In the simulation, we used M = 50 



in all the TSP cases that we tested. We performed 150 
Monte Carlo steps in each layer at a temperature that is 
near the optimal temperature in each of the cases which 
we give in Table I. We have tested five cases from the 
TSP databank [12]: dl98, lin318, pcb442, rat 783 and 
A1577. These cases are selected from the set of sample 
problem instances as listed in [1]. We have performed 100 
trial runs in each case and computed the average value 
of the best solutions obtained in each of the 100 trial 
runs. We also kept record of the best solution among 
the 100 trial runs and list them in Table I. We used 6 
layers for the first three cases and 8 layers for the last 
two cases respectively. To make our discussion simple, 
we only consider 2-opt moves in our algorithm. Other 
moves can be considered but is not our main emphasis 
here. 

Table I. Tests on 5 cases from TSPLIB. T is the temperature used 
for the Monte Carlo simulation, RQ([1]=% above the optimal solu- 
tion) is the relative quality of the average length of the best paths 
obtained for the 100 trial runs with respect to the optimal solutions 
given in TSPLIB, Best is the shortest path length obtained in the 
100 trial runs using the present algorithm and Op is the optimal 
solution given in TSPLIB. 



Case 


T 


RQ 


Best 


Op 


dl98 


10.0 


0.0155% 


15780 


15780 


lin318 


37.5 


0.30% 


42029 


42029 


pcb442 


27.5 


0.67% 


50798 


50778 


rat783 


4.0 


0.68% 


8819 


8806 


A1577 


2.75 


0.88% 


22273 


22249 



As can be seen in Table I, the algorithm can locate 
the optimal paths in the two cases with fewest cities, 
i.e. dl98 and lin318 while it can obtain near-optimal 
solutions for the other three cases within the number of 
Monte Carlo steps we preset. In general, the more the 
individual Monte Carlo runs, i.e. the larger the M, and 
the more the Monte Carlo steps and layers one uses, the 
better will be the solution obtained. For example, in the 
case of dl98, 23 out of 100 trial runs found the optimal 
solution when we perform 150 Monte Carlo steps in each 
layer. When we increased the Monte Carlo steps in each 
layer to 200, the average value of the best runs would 
now be only 0.00906% above the optimal solution with 
41 out of 100 trial runs found the optimal solution. As a 
comparison, we have performed a simulation on the dl98 
case using the conventional Monte Carlo method, for 100 
trial runs with 45,000 Monte Carlo steps per run. The 
number of Monte Carlo steps here is equivalent to the 
number of M (50 here) times the number of layers (6 
in the case of dl98) times the 150 Monte Carlo steps in 
each layer. The average value for the best solution of the 
100 trial runs is 1.07% above the optimal value and the 
best solution is 15843. A further increase of Monte Carlo 
steps to 10 5 improves the average value slightly to 0.90% 
above the optimal value. This shows the superiority of 



the present algorithm over the conventional Monte Carlo 
method. 

It is perhaps worth mentioning that locating the op- 
timal or near-optimal temperature should help one per- 
form simulations in the case of random TSP in a two- 
dimensional unit square or something similar. As an ex- 
ample, we have obtained an average optimal tempera- 
ture around 0.02 over 10 configurations of the random 
TSP with 100 cities, when 100 Monte Carlo steps are 
performed in each layer of a 6 layer simulation. The opti- 
mal temperature range is however affected by the number 
of Monte Carlo steps and the moves(we use 2opt moves 
here) . A rescaling of the optimal temperature range with 
respect to the system size, i.e. the number of cities, is 
also needed. 

(2) The Spin Glass Problem (SGP) 

Spin glasses have been a subject in statistical physics 
that is under intensive study. The simplest spin glass 
model is the Edwards-Anderson model and its Hamil- 
tonian has been given in the above. There are exact 
methods available [13] for such spin glass models in two- 
dimensions but are not effective in three-dimensions. Es- 
timates of the ground state energy are available for the 
three dimensional case [14,15] in the literature for lattice 
sizes up to L = 12 and could be used for comparison. 
Our task here is to use the present algorithm to esti- 
mate the ground state energy of the three-dimensional 
Anderson-Edwards model within a reasonable amount of 
CPU time. Employing our algorithm to this model is 
straightforward. To begin with, we have a set of M ran- 
domly generated spin configuration for an L x L x L 
lattice with randomly assigned ± J couplings. We again 
let each of the configurations evolve for a preset Monte 
Carlo steps and record the M best solutions from the M 
individual runs. These M best solutions are then used 
to construct the guided distribution function. Similar 
to the TSP, we freeze a sub-spin configuration pattern 
on this three-dimensional lattice if it exceeds a certain 
preset probability. In this way, we reduce the spin de- 
grees of freedom after each layer of simulation. We again 
randomly generate M spin configurations at the begin- 
ning of each new layer. There is one thing one needs 
to keep in mind in constructing the guided distribution 
function of the spin glass problem. When all the spins 
on the three-dimensional lattice changes sign, the energy 
of the system does not change. Therefore, one needs to 
fix one of the spins to take on the same value for all the 
M best solutions when one constructs the guided distri- 
bution function. For example, we can set the spin value 
of the (0,0,0) site of the L x L x L lattice to be 1. If a 
particular best solution has a value -1 for this site, then 
all the spins have to flip sign in this best solution before 
one records its spin configuration. 

In our simulation, we have used M = 15. We have 
carried out simulations for L ranging from 4 to 12 and 
the result is given in Table II. Since our purpose here 
is to obtain the ground state energy of the model, we 



will instead perform several runs for each spin configura- 
tions and pick the lowest energy as its ground state. The 
number of layers and the number of Monte Carlo steps 
in each layer that we used will also depend on L. For ex- 
ample, for L — 6, we used 2 layers with 150 Monte Carlo 
steps in each layer and 5 trial runs for each randomly as- 
signed coupling configuration we generated. In the case 
of L = 8 sites, we used 6 layers with 200 Monte Carlo 
steps in each layer and 10 trial runs for each randomly 
assigned coupling configuration before we take the best 
solution as the ground state energy for that particular 
spin configuration. The number of Monte Carlo steps 
and layers for different L are given in Table II. 

In the simulation of the spin glass model, we have ob- 
served that unlike the case of the TSP, the temperature 
range in the Boltzmann factor that is most efficient to 
locate the ground state or near ground state energy is 
between T = J and T = 1.1J for all the L that we tested 
and we have used T — 1.05 J for all our simulation of the 
spin glass model. That the optimal temperature range is 
not sensitive to the lattice size here is probably related 
to the fact that one only considers the nearest neighbor 
interaction in the model. 



speed up of optimum search processes, it proves to be a 
powerful new algorithm. It takes about 24 seconds and 9 
seconds of a Pentium II 450 MHz CPU time to perform 
a trial run of the dl98 of the TSP and L = 8 of the SGP. 
Tests of our algorithm in the TSP and SGP are promis- 
ing. Incorporating heuristic methods into this algorithm 
will certainly improve the efficiency and power of its op- 
timum search. This algorithm is very general and can 
be applied to many optimization problems including dis- 
crete optimization problems we tested here and optimiza- 
tion problems with continuous parameters. The present 
algorithm has been applied to the Lennard- Jones micro- 
cluster problem and can reproduce all the known ground 
state energy for the system sizes that the algorithm has so 
far tested [17] . We believe that this algorithm should be a 
simple and general enough algorithm that can obtain op- 
timal or near-optimal solutions within a reasonable CPU 
time for optimization problems. 



Table II. Tests on L X L X L SG instances, m is the number of 
randomly assigned ±J (J = 1 here) cases, MC is the number of 
Monte Carlo steps in each layer and E is the average ground state 
energy of the M spin configurations. We have set J = 1 in our 
simulation. We also include here the results from Ref[14,15] for 
comparison. 



/;? 



MC 



E 



Ref[14] Ref[15] 



4 200000 50 

6 50000 150 

8 4000 200 

10 400 400 

12 40 600 



-1.73750(13) 
-1.77128(12) 
-1.77910(28) 
-1.7831(8) 
-1.7847(14) 



-1.73749(8) -1.7370(9) 
-1.77130(12) -1.7723(7) 
-1.77991(22) -1.7802(5) 
-1.78339(27) -1.7840(4) 
1.78407(121) -1.7851(4) 



As shown in Table II, the ground state energy that we 
obtained by using our algorithm is comparable with the 
present estimate of the ground state energy of the model 
by using other methods. If more layers, Monte Carlo 
steps and trial runs are used, the result will be improved. 
In this paper, we have introduced a new algorithm by 
incorporating a guided distribution function into the con- 
ventional Monte Carlo method. This new algorithm is 
both biologically and physically motivated. By evolving 
a population of M solutions and record their best individ- 
ual solutions in each layer, we are relying on the cooper- 
ative effort of the population, similar to many evolution- 
ary algorithms such as genetic algorithms [16]. Reducing 
the computational phase space by eliminating irrelevant 
degrees of freedom of the problem as we continue our 
layer by layer simulation is similar to the renormalization 
group equation method in statistical physics. Together 
with the use of an optimal temperature in Monte Carlo 
simulation for optimization problems, which results in a 
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